Stability of Finite-Difference Models Containing Two Boundaries or Interfaces

It is known that the stability of finite-difference models of hyperbolic initial-boundaq value problems is connected with the propagation and reflection of parasitic waves. Here the waves point of view is applied to models containing two boundaries or interfaces, where repeated reflection of trapped wave packets is a potential new source of instability. Our analysis accounts for various known instability phenomena in a unified way and leads to several new results, three of which are as follows. (1) Dissipativity does not ensure stability when three or more formulas are concatenated at a boundary or internal interface. (2) Algebraic "GKS instabilities" can be converted by a second boundary to exponential instabilities only when an infinite numerical reflection coefficient is present. (3) "GKS-stability" and "P-stability" can be established in certain problems by showing that the numerical reflection coefficient matrices have norm less than one. 0. Introduction. Hyperbolic systems of partial differential equations admit solutions which behave locally like waves moving along characteristics. When such a system is modeled numerically by finite differences or finite elements, the result is a dispersive medium that may admit additional parasitic wave modes with wave-lengths on the scale of the discretization. Energy associated with these parasitic waves travels at a group velocity that is unrelated to the characteristics of the original system [25], [30]. However, the behavior of such waves has a decisive effect on stability. For finite-difference models of linear hyperbolic problems with a single spatial boundary, a stability theory was developed around 1970 by Kreiss, Osher, and others [lo], [18]. In earlier papers we have shown that this theory can be naturally stated in terms of dispersive wave propagation [26], [27]. To summarize: if a boundary with homogeneous boundary conditions can emit a radiated wave in the absence of any incident waves, i.e., a wave with group velocity pointing into the interior of the domain, then it is unstable. If it has an infinite reflection coefficient for waves at some frequency, a stronger condition, then it is more severely unstable. This paper applies wave propagation ideas to investigate stability for one-dimensional linear finite-difference models with two or more boundaries or internal interfaces. The most basic example of such a model is a discrete approximation to an equation whose spatial domain is an interval such as [O,l]. Another example is a model of a problem featuring discontinuous coefficients, e.g., wave propagation in a discontinuously stratified medium [4]. A third is a numerical scheme employing local Received October 26, 1984. 1980 Mathematics Subject Classification. Primary 65M10. *Research was supported in part by the U. S. Department of Energy under contract DE-AC02-76ER03077-V and in part by the National Aeronautics and Space Administration under NASA Contract No. NAS1-17070 while the author was in residence at the Institute for Computer Applications in Science and Engineering, NASA Langley Research Center, Hampton, VA 23665. C1985 American Mathematical Society 0025-5718/85 $1.00 + $.25 per page 280 LLOYD N. TREFETHEN mesh refinement to improve accuracy, in which various interfaces between fine and coarse grids will be present [2], [5]. A fourth is any model with a composite numerical boundary or interface, such as a fourth-order difference model on [0, m ) that has a five-point stencil, and whch therefore requires one numerical boundary condition at j = 0 and another at j = 1 [15]. Such a model can be viewed as containing two interfaces separated by a single grid point, and we will show that this view may be useful for stability analysis. Any multi-interface model potentially admits trapped numerical waves that reflect back and forth repeatedly from one interface to another. If the reflections cause amplification, this can lead to unbounded growth of numerical solutions. The factors that control this are: magnitude of the reflection coefficients, which is the source of growth; dissipation of waves as they travel between interfaces, which is a source of attenuation; and travel time between interfaces, whch determines how frequently any reflection circuit that causes growth is repeated. All of the arguments of this paper consist of worlung out consequences of various combinations of these factors that may be of practical interest. In particular, we investigate two kinds of stability. First, "stability" or "LaxRichtmyer stability" refers to the usual Lax-Richtmyer definition for time-dependent finite-difference models, or to variants thereof such as "GKS-stability" (Definition 3.3 in the well-known paper by Gustafsson, Kreiss, and Sundstrom [lo]). A difference model that is stable in this sense may admit solutions that grow with time, provided that the growth does not get worse as the mesh is refined. This is what is needed to ensure convergence as the mesh size approaches zero to the correct solution of the time-dependent differential equation, for each fixed time t . On the other hand, to be "P-stable" [I], a model must admit no growing solutions at all. (See Section 3 for the precise definition. Such a model is also sometimes called "time-stable" [29].) Although the theory here is not as well developed, such a condition is needed if a time-dependent difference model is to be used to obtain steady-state solutions, as is common in computational aerodynamics. As a rule of thumb, we will show that P-instability is very often associated with reflection coefficients greater than 1 in magnitude, and Lax-Richtmyer instability with reflection coefficients that are infinite. Section 1 reviews stability theory for one-boundary difference models (Proposition 1). Section 2 investigates interfaces separated by a fixed number of grid points A j as the mesh is refined, as in the fourth-order boundary condition mentioned above. Here the travel times go to zero with the mesh spacing, with the effect that finite reflection coefficients greater than 1 can cause catastrophic unstable growth (Propositions 3, 4, 4'). Conversely, reflection coefficients smaller than 1ensure stability (Proposition 5). Section 3 considers interfaces separated by a fixed distance Ax as the mesh is refined, as in the problem on [O,l] mentioned above. Here, the travel times are independent of mesh spacing, so large finite reflection coefficients can cause P-instability (Proposition 7), but not Lax-Richtmyer instability (Proposition 6). In this context multiple reflections may convert the weak instability of a single interface to a catastrophic instability (Proposition 8), but only if the unstable interface is of the sort with an infinite reflection coefficient (Proposition 9). Once again, reflection coefficients smaller than 1in norm ensure stability (Proposition lo), 281 STABILITY OF FINITE-DIFFERENCE MODELS and if the model is dissipative, it suffices to look at the reflection coefficients for the partial differential equation rather than its numerical approximation (Proposition 11). For convenience of reference, here is a list of the explicit examples presented here to illustrate various points. The symbol A indicates the modulus of a reflection coefficient, and S, the solution norm at time step n. These quantities will be made more precise later on. 1. Algebraically unstable one-boundary model (one boundary, A = co, S, const n). 2. Exponentially unstable concatenation of three stable dissipative formulas (fixed-Aj, A > 1,Snconstn). 3. P-stability guaranteed by reflection coefficients less than 1 (fixed-Aj or Ax, A G 1, Snconst). 4. P-instability caused by reflection coefficients greater than 1 (fixed-Ax, A > 1,S, const'). 5 . Exponential instability caused by interaction of two algebraically unstable boundaries (fixed-Ax, A = co,Sn(A j)co"St'). The reader may be disappointed at the artificiality of some of these examples, especially (2) and (3), and he may wonder how helpful wave reflection ideas can be in practice for the design of difference schemes. A full answer to this question will have to await further experience. Nevertheless, there is no doubt that the instability mechanisms described here are real and deserve to be understood. At present, virtually no difference models containing multiple interfaces have been shown to be stable. Perhaps the ideas here, such as Proposition 5 , can help bring about a change in this situation. Much of the material in thls paper can be found in Section 6 of the author's Ph. D. dissertation [24]. For some numerical illustrations, see [28]. For a different analysis of stability for two-boundary problems that is closely related to the present one, see the report [8] by Giles and Thompkins, which is mainly concerned with P-stability. Giles and Thompluns consider parasitic wave propagation for models with variable as well as constant coefficients. The phenomenon of instability caused by trapped wave packets can also occur in two-dimensional problems when the domain contains a corner. Osher has given examples of hyperbolic systems (not difference models) in corners that are ill-posed because of trapped waves [19], while Sarason and Smoller have shown that for a 2 x 2 strictly hyperbolic system such as the second-order wave equation, this cannot happen [21]. But trapped numerical waves may render a finite-difference model of even the latter sort unstable. The principles involved are precisely those of this paper, but we will discuss corners elsewhere. The reader interested in getting to the main ideas quickly may find it possible to turn directly to Section 2. 1. Review of Wave Propagation and Stability for One-Boundary Difference Models. Consider a linear first-order hyperbolic system of partial differential equations